home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / anova.pro < prev    next >
Text File  |  1997-07-08  |  14KB  |  465 lines
  1. ; $Id: anova.pro,v 1.3 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ;  Copyright (c) 1991-1997, Research Systems Inc.  All rights
  4. ;  reserved. Unauthorized reproduction prohibited.
  5.  
  6.  
  7.     Pro ContrastH,T,M,SS, MS, Names, n, rep, DFE, unit
  8.  
  9. ;ContrastH evaluates the significance of the
  10. ;contrasts contained in the 1-or 2- dim array T and
  11. ;inserts the results into the anova table under
  12. ;construction. T(i,j)= coefficient of ith mean in
  13. ;jth contrast. M/S= means to be used in contrasts,
  14. ;MS = mean square error.  Names= treatment,block,or
  15. ;interaction names n=0,1,2,3 depending on whether
  16. ;contrast is of 1-way treatments 2-way with/without
  17. ;interactions.Rep = # of replications
  18.  
  19.  
  20.  S = SS
  21.  SM=size(M)            ;Size Info for Means
  22.  if(SM(0) EQ 1) THEN  S1= SM(1) else S1= SM(1)*SM(2)
  23.  
  24.  
  25.  ST=size(T)                     ; Size info for Coefficients
  26.  C=ST(1)
  27.  
  28.  OK = testcontrast(T, unit)
  29.  if OK EQ 0 THEN BEGIN
  30.    printf,unit,'anova ---contrast array contains equations'
  31.    printf,unit,'whose coefficients do not sum to 0'
  32.    return
  33.  ENDIF
  34.  if OK EQ -1 then return
  35.  
  36.  if(ST(0) GT 1) Then Con=ST(2) else Con=1
  37.  
  38.  SN=Size(Names)
  39.  if(C NE S1 ) THEN BEGIN        ; Check compatibility
  40.    printf,unit, 'anova- Contrast has wrong size'
  41.    Return
  42.  ENDIF
  43.  
  44.   if(SN(0) EQ  0) THEN Name="Contrast" +  $
  45.                             StrTrim( INDGEN(Con),2) $
  46.   else if   (SN(1) LT Con) THEN BEGIN
  47.    Nm=INDGEN(Con)
  48.    Name=[Name,"Contrast" +StrTrim(Nm(SN:Con-1))]
  49.    ENDIF else Name=Names
  50.  
  51.  
  52.  if(ST(0) EQ 1)  then BEGIN 
  53.                             ;comput SS for each contrast
  54.      P = Total(T*M)
  55.      if(n EQ 3) THEN S=1
  56.      P = P^2/Total(T*T)
  57.      P = P/(S*rep)               
  58.  
  59.     printf,unit,           $
  60.     Format=                $
  61. '(5X,A17,1X,G13.4,3X,I5,3X,G11.4,1x,G11.4, 3X, G11.4)',$
  62.              Name(0),P,1,P,P/MS, 1-F_Test1(P/MS,1,DFE) 
  63.  
  64.  ENDIF else BEGIN
  65.               R1=ST(2)
  66.               P=Fltarr(R1)
  67.               A=reform(M/S,C)
  68.               if(n EQ 3) THEN S=1
  69.               P = A # T
  70.              T2 = replicate(1.0,C) #  T^2 
  71.              P = S * rep * P^2 / T2
  72.               for i =0L,R1-1 DO BEGIN
  73.                 printf,unit,                      $
  74.                 Format=     $
  75.                 '(5X,A17,1X,G13.4,3X,I5,3X,G11.4,1X,G11.4,3X,G11.4)', $
  76.                            Name(i), P(i),1, P(i),P(i)/MS,    $
  77.                                  1-F_Test1(P(i)/MS,1,DFE)
  78.               ENDFOR
  79.        ENDELSE
  80.  
  81.   RETURN
  82.   END                      ;ContrastH 
  83.  
  84.                                  
  85.                                  
  86.                                  
  87.                                                                      
  88.                                
  89.  
  90.  PRO  Anova, X1, FCTEst = FCTEst, FRTest = FRTest,    $
  91.                 FRCTest = FRCTest, DFE = DFE, DFC = DFC, $
  92.                  DFR = DFR, DFFRC = DFRC,   $
  93.             One_Way = One, Unequal_One_Way = Unequal,$
  94.              Two_Way = Two, Interactions_Two_Way  = $
  95.             Interactions, Missing=M, List_Name=LN, $
  96.             TContrast=TC, BContrast=BC, IContrast=IC, TName=TN,$
  97.             BName=BN,TCName=TCN, BCName=BCN, ICName=ICN,  $
  98.             No_Printout = NP
  99. ;+
  100. ;
  101. ; NAME:
  102. ;    ANOVA
  103. ;
  104. ; PURPOSE:
  105. ;    Perform one-way analysis of variance with
  106. ;    equal and unequal sample size, two-way analysis
  107. ;    of variance with and without interactions.
  108. ;
  109. ; CATEGORY:
  110. ;    Statistics.
  111. ;
  112. ; CALLING SEQUENCE:
  113. ;    ANOVA, X
  114. ;
  115. ; INPUTS: 
  116. ;    X:    An array of experimental data. For one-way and two-way 
  117. ;        analysis of variance without interactions, X must be
  118. ;        dimensioned (Treatment#,B),  where B is the maximum
  119. ;        sample size (one-way anova) or number of blocks
  120. ;        (two-way analysis).  With interactions, X is dimensioned 
  121. ;        (Treatment#,I,B), where I is the number of iterates in each
  122. ;        cell.
  123. ;             
  124. ; OUTPUT:
  125. ;    Anova table displaying Sum of Squares, Mean Squares,
  126. ;    and F Ratio for sources of variation.
  127. ;
  128. ; KEYWORDS:
  129. ;    FCTest = if present and set, returns the value of F for
  130. ;        treatment or column variation.
  131. ;    FRTest    if present and set, returns value of F for row
  132. ;        or block variation.
  133. ;    FRCTest    if present and set, returns value of F for column
  134. ;        variation.
  135. ;    DFE    if present and set, returns denominator degrees of
  136. ;        freedom for FCTest,FRTest,FRCTest.
  137. ;    DFC    if present and set, returns numerator degrees of
  138. ;        for FCTest 
  139. ;    DFR    if present and set, returns numerator degrees of for
  140. ;        FRTest. 
  141. ;    DFRC    if present and set, returns numerator degrees of
  142. ;        for FRCTest. 
  143. ;
  144. ;    Missing    missing data value. If undefined,
  145. ;        assume no missing data.If unequal
  146. ;        sample sizes, set M to place holding
  147. ;        value.
  148. ;
  149. ;    List_Name:    name of output file. Default is to the
  150. ;        screen.
  151. ;
  152. ;  Type of design structure. Options a
  153. ;    ONE_WAY         =  if set, perform 1-way anova.
  154. ;    Unequal_ONE_WAY =  if set, perform 1-way ANOVA with
  155. ;                             unequal sample sizes.
  156. ;    TWO_WAY = if set, perform a 2-way ANOVA.
  157. ;    Interactions_TWO_WAY = if set, perform a
  158. ;                                2-way ANOVA with
  159. ;                                interactions.
  160. ;  One and Only one of the options above must be set.
  161. ;
  162. ;    TContrast:    1- or 2- dimensional array 
  163. ;        of treatment contrasts.
  164. ;        Each row of TC represents a contrast.
  165. ;        TC(i,j) = the coefficient of the mean of
  166. ;        the ith treatment in the jth contrast. 
  167. ;
  168. ;    BContrast:    1- or 2- dimensional array of block
  169. ;        contrasts. Each row of BC represents
  170. ;        a contrast. BC(i,j) = the coefficient
  171. ;        of the mean of the ith block in the jth
  172. ;        contrast.                 
  173. ;
  174. ;    IContrast:    1- or 2- dimensional array of interaction contrasts.  Each 
  175. ;        row of TC represents a contrast.  IC(i,j) = the coefficient of
  176. ;        the mean of the ith treatment in the jth contrast. 
  177. ;
  178. ;    TName:    name to be used in the output for treatment type.
  179. ;
  180. ;    BName:    name to be used in the output for block type.
  181. ;    TCName:    vector of names to be used in the output to identify 
  182. ;        treatment contrasts.
  183. ;    BCName:    vector of names to be used in the output to identify block
  184. ;        contrasts.
  185. ;    ICName:    vector of names to be used in the output to identify
  186. ;        interaction contrasts.
  187. ;  No_Printout:    flag, when set, to suppress printing of output.
  188. ;
  189. ; RESTRICTIONS:
  190. ;    NONE.
  191. ;
  192. ; SIDE EFFECTS:
  193. ;    None.
  194. ;
  195. ; PROCEDURE:
  196. ;    Calculation of standard formulas for sum of squares, mean squares and 
  197. ;    F ratios.        
  198. ;  
  199. ;
  200. ;-
  201.  
  202.  On_Error,2
  203.  if N_elements(X1) NE 0 THEN X= double(X1)
  204.  SX=Size(X)        
  205.  C=SX(1)      ; Compute number of columns       
  206.  D= 1         ; Default # of replications
  207.  
  208.  if (N_ELements(LN) EQ 0) THEN  Unit = -1      $
  209.  ELSE openw,unit,/Get,LN
  210.  
  211.  ONE_WAY = 0 & ONE_WAY_UNEQUAL=1 & TWO_WAY = 2
  212.  TWO_WAY_INTERACTIONS =3
  213.  
  214.  keyset = [keyword_set(One), keyword_set(Unequal), $
  215.            keyword_set(Two),keyword_set(Interactions)]
  216.   T  = where(keyset,nkey)
  217.  
  218.  if nkey NE 1 THEN BEGIN
  219.     printf,unit,    $
  220. 'anov - must specify one and only one type of design structure'
  221.     goto, DONE
  222.     ENDIF     
  223.  
  224.  T = T(0)
  225.  if((T EQ TWO_WAY_INTERACTIONS) AND  (SX(0) NE 3)) THEN BEGIN
  226.     printf,unit, 'Anova- wrong number of dimensions' 
  227.     goto,DONE
  228.     ENDIF
  229.  
  230.  if T NE TWO_WAY_INTERACTIONS THEN  R=SX(2) ELSE BEGIN
  231.    R=SX(3) 
  232.    D=SX(2)    
  233.  ENDELSE
  234.                                        ;Compute number of rows
  235.  
  236.  
  237.  if(C LE 1 or R LE 1 ) THEN BEGIN 
  238.    printf,unit,   $
  239.          'ANova- all dimensions must be greater than one' 
  240.    goto,Done
  241.    ENDIF
  242.  
  243.  if(N_Elements(TN) NE 0) THEN TR=TN ELSE TR='Column'
  244.                                     ;treatment names 
  245.  
  246.   if(N_Elements(BN) NE 0) THEN Bl=BN else BL='Row'
  247.                                 ;block names
  248.  
  249.  
  250.  
  251.  
  252.  
  253.  if( T EQ ONE_WAY_UNEQUAL and N_Elements(M) ne 0) THEN BEGIN
  254.  
  255.     
  256.      Pairwise,X,M,YR,Y,notgood,good           ; replace missing data by 0 
  257.  
  258.      if(R LT  Max(Y)) THEN BEGIN   ; still enough rows ? 
  259.    printf,unit,      $
  260. 'Anova Row number of Data array too small for population sizes'  
  261.      goto,DONE
  262.    ENDIF
  263.  
  264.      if (Min(Y) LE 1) THEN BEGIN
  265.        printf,unit,       $
  266.              ' anova- each sample size must be greater than 1'
  267.        goto,DONE
  268.       ENDIF
  269.      
  270.      N = Total(Y)
  271.  
  272.  ENDIF ELSE N = R * C *D  
  273.  
  274.  
  275.  
  276.  MEAN = Total(X)/N           ; compute mean
  277.  
  278.          ;compute total sum of squares
  279.  if T EQ ONE_WAY_UNEQUAL and N_ELEments(Missing) ne 0 THEN $
  280.     SST = Total((X(good) - mean)^2)   $
  281.  else SST= Total((X-mean)^2)    
  282.  
  283.  
  284.     ; Pre-process three dimensional array by converting toa
  285.     ; a two - dimensional array whose entries are sums of
  286.     ; repetitions 
  287.  
  288.  if(T EQ TWO_WAY_INTERACTIONS) THEN BEGIN    
  289.      X= dblarr(C,R)
  290.        
  291.     for i = 0L,R-1 Do for j=0L,C-1 DO          $
  292.      X(j,i) = Total(X1(j,*,i)) 
  293.                      
  294.   ENDIF 
  295.  
  296.  
  297.  
  298.  
  299.  ColTotal=X # replicate(1.,R)   ; compute vector of column totals
  300.  
  301.  
  302.  
  303.  DFR= long(R-1)                     ; row degrees of freedom
  304.  DFC= long(C-1)                     ; column degrees of freedom
  305.  
  306.  case 1 of
  307.               
  308. T EQ ONE_WAY_UNEQUAL OR T EQ ONE_WAY  :       $
  309.       Begin              
  310.             ;compute treatment sum of squares
  311.  
  312.             if T EQ ONE_WAY_UNEQUAL and N_elements(M) ne 0 THEN      $
  313.                SStr = Total(y*(Coltotal/Y - mean)^2) $
  314.             else SStr = Total(R*(Coltotal/R - mean)^2) 
  315.  
  316.             SSE = SST - SStr               ;compute error sum of squares
  317.             DFT=N-1                     
  318.             DFE=N-c                        ; error degrees of
  319.                                            ; freedom
  320.             
  321.             MSSE=SSE/DFE                   ; error mean square
  322.             if(MSSE EQ 0) THEN BEGIN
  323.                printf,unit,                $
  324.     'anova- must stop since mean square error =0'
  325.                goto,DONE
  326.             END
  327.             END
  328.  
  329.  
  330.  
  331. T EQ TWO_WAY:Begin 
  332.             RowTotal = replicate(1,C) # X     
  333.                               ;computations for two-way anova
  334.             SSR=  C*Total((RowTotal/C - mean)^2)
  335.             SSTr= R *Total((ColTotal/R - mean)^2)
  336.             SSe=SST-SSTr-SSR
  337.             MSSR=SSR/DFR                      ;row mean square
  338.             DFE=(r-1)*(c-1)
  339.             MSSE=SSE/((r-1)*(c-1))
  340.             if(MSSE EQ 0) THEN BEGIN
  341.             printf,unit,                      $
  342.                  'anova- must stop,since mean square error =0'
  343.             goto,DONE
  344.             END
  345.             FRTest=MSSR/MSSE                 ;F value for rows
  346.             DFT=R*C-1
  347.             END
  348.  
  349.  T EQ TWO_WAY_INTERACTIONS:BEGIN
  350.             RowTotal = Replicate(1./D,C)#X
  351.                     ;computations for two-way with interactions
  352.             ColTotal = ColTotal/D
  353.             MR = RowTotal/(C)    ; Row means
  354.             MC = ColTotal/(R)    ; Column means
  355.             MI = X/D;              ; Interaction means
  356.  
  357.             MR1 = Replicate(1.,C) # MR
  358.             MC1 = MC # Replicate(1.,R)
  359.  
  360.                  ; Treatment sum of squares         
  361.             SSTr= D * R * Total((MC - mean)^2)
  362.  
  363.                  ; Block sum of squares 
  364.             SSR = D * C * Total((MR - mean)^2)
  365.  
  366.                 ; Interaction sum of squares
  367.            SSRC = D * Total((MI - MR1 - MC1 + mean)^2)
  368.             
  369.  
  370.             SSE=SST-SSR-SSTR-SSRC
  371.             MSSR=SSR/DFR
  372.             DFRC=(r-1)*(C-1)
  373.             MSSRC=SSRC/DFRC
  374.             DFE=r*c*(D-1)
  375.             MSSE= SSE/DFE
  376.             if(MSSE EQ 0) THEN BEGIN
  377.             printf,unit,        $
  378.                  'anova-must stop,since mean square error =0'
  379.             goto,DONE
  380.             END
  381.             FRTest=MSSR/MSSE
  382.             FRCTest=MSSRC/MSSE
  383.             DFT=r*c*d-1            
  384.             End
  385.  
  386.             
  387.  
  388.  
  389.  
  390.   Else: BEGIN
  391.         printf,unit,     $
  392.     "anova - keyword Structure has unknown value"
  393.         goto,done
  394.         end
  395.  ENDCASE
  396. ;
  397. ;
  398.       MSSC=SSTr/(c-1) 
  399.                     ; Column mean square for all types of anova
  400.       FCTest=MSSC/MSSE           
  401.                      ;Column F value for all types of anova
  402.   if (N_Elements(NP) EQ 0) THEN BEGIN
  403.     
  404.  
  405.  
  406.      
  407.   
  408.       printf,unit,'          Source        SUM OF SQUARES    DF     MEAN SQUARE       F       p'
  409.       printf,unit,'          **************************************************************************'
  410.  
  411.  
  412.     if (T EQ TWO_WAY) OR (T EQ TWO_WAY_INTERACTIONS) THEN BEGIN
  413.       printf,unit,       $
  414.        Format=    $
  415.   '(A14,7X,G15.7,3X,I5,3X,G11.4,1X,G11.4,3X,G11.4)',$
  416.                 Bl,SSR,DFR,MSSR,FRTEst, 1-F_Test1(FRTEST,DFR,DFE)
  417.  
  418.       if (N_Elements(BC) NE 0) THEN                 $
  419.          ContrastH,BC,RowTotal,C,MSSe,BCN,2,D,DFE,unit 
  420.     ENDIF
  421.  
  422.       printf,unit," "
  423.       printf,unit,                                     $
  424.  
  425.      Format='(A14,7X,G15.7,3X,I5,3X,G11.4,1X,G11.4,3X,G11.4)',$
  426.              Tr,SSTR,DFC,MSSC,FCTEst, 1-F_Test1(FCTEST,DFC,DFE)
  427.  
  428.     if(N_Elements(TC) NE 0) THEN if(T NE ONE_WAY_UNEQUAL)   $
  429.       THEN ContrastH,TC,ColTOtal,R, MssE,TCN,0,D,DFE,unit
  430.       
  431.        printf,unit," "
  432.  
  433.  
  434.       if T EQ TWO_WAY_INTERACTIONS THEN BEGIN 
  435.       printf,unit,                    $
  436.       Format=      $
  437.       '(A14,7x,G15.7 ,3x,I5,3x,G11.4,1x,G11.4,3x,G11.4)', $
  438.       'Interaction',SSRC,DFRC,MSSRC,FRCTest,         $
  439.                                    1-F_Test1(FRCTEST,DFRC,DFE)
  440.       if(N_Elements(IC) NE 0) THEN    $
  441.         ContrastH,IC,X,D,MSSe,ICN,3,D, $
  442.                                               DFE,unit
  443.        ENDIF
  444.       printf,unit, " "
  445.       printf,unit,         $
  446.        Format='(A14,7X,G15.7,3X,I5,3X,G11.4)',         $
  447.                                'Error',SSe,DFe,MSSe
  448.       printf,unit, " "
  449.       printf,unit,    $
  450.           Format='(A14,7x,G15.7,3X,I5)', 'Total',SST,DFT
  451.      ENDIF
  452.  
  453.  DONE:
  454.    if ( unit NE -1) THEN Free_Lun,unit
  455.    Return
  456.  END
  457.             
  458.  
  459.  
  460.  
  461.             
  462.  
  463.  
  464.  
  465.